Communities of endophytic fungi in a Puerto Rican rainforest vary along a gradient of disturbance due to Hurricane Maria

Abstract Increases in the frequency and intensity of hurricanes influence the structure, function, and resilience of Caribbean forests. Trees in such forests harbor diverse fungal endophytes within leaves and roots. Fungal endophytes often are important for plant health and stress responses, but how their communities are impacted by hurricanes is not well known. We measured forest disturbance in Carite State Forest in Puerto Rico ca. 16 months after the passage of Hurricane Maria, a Category 4 storm. In three sites, each comprising three plots representing a local gradient of hurricane disturbance, we evaluated soil chemistry and used culture‐free analyses to measure richness, phylogenetic diversity, and composition of endophyte communities in leaves and roots. We found that endophyte richness did not vary significantly among plant families or as a function of soil chemistry. Instead, leaf endophytes peaked in richness and decreased in phylogenetic diversity at intermediate levels of disturbance. Root endophytes did not show such variation, but both leaf‐ and root endophyte communities differed in species composition as a function of disturbance across the forest. Locations with less disturbance typically hosted distinctive assemblages of foliar endophytes, whereas more disturbed locations had more regionally homogeneous endophyte communities. Together, our results show that changes in endophyte richness and phylogenetic diversity can be detected in aboveground tissues more than a year after major storms. In turn, pervasive shifts in endophyte community composition both aboveground and belowground suggest a subtle and lasting effect of hurricanes that merits further study, potentially contributing to the promotion of spatially heterogeneous endophyte assemblages at a landscape scale in these diverse island forests.

On September 20, 2017, Hurricane Maria struck Puerto Rico as a Category 4 storm with sustained winds up to 249.5 km/h and heavy rainfall (López-Marrero et al., 2019;Pasch et al., 2018;Wachnicka et al., 2020). As the strongest hurricane in Puerto Rico in nearly 90 years, Maria's wind gusts, rainfall, and landslides damaged mature forests in many parts of the island, with an estimated 23-31 million trees severely damaged or killed (Feng et al., 2018; see also Uriarte et al., 2019). Accordingly, comparisons of pre-hurricane (2015 and 2016) and post-hurricane surveys revealed reductions in vegetation cover after Hurricane Maria's passage (Hosannah et al., 2020;Hu & Smith, 2018;Van Beusekom et al., 2018).
For these reasons, endophytes are increasingly considered important for forest health (Arnold et al., 2003;Bamisile et al., 2018;Christian et al., 2019;Griffin & Carson, 2018;Rodriguez et al., 2009;Yamaji et al., 2016), and understanding their responses to disturbances such as those caused by hurricanes is important for understanding the dynamics of forest regeneration (see Busby et al., 2022).
Effects of hurricanes on soilborne microbial communities in forests have been studied previously (e.g., Cantrell et al., 2014;Eaton et al., 2020), providing some insight into the ways in which fungi diversos hongos endófitos dentro de las hojas y raíces. Los hongos endófitos pueden ser importantes para la salud de las plantas y las respuestas a el estrés, pero no se conoce bien cómo los huracanes afectan a sus comunidades. Medimos la perturbación forestal en el Bosque Estatal de Carite en Puerto Rico a 16 meses después del paso del Huracán María (huracán categoría 4). En tres sitios, cada uno de los cuales incluyeron tres parcelas que representan un gradiente local de perturbación debido al huracán, evaluamos la química del suelo y utilizamos análisis libres de cultivos para medir la riqueza, la diversidad filogenética y la composición de las comunidades de endófitos en las hojas y raíces de las plantas. Encontramos que la riqueza de endófitos no varió significativamente entre las familias de plantas o en función de la química del suelo.
Los endófitos de hojas alcanzaron su punto máximo en riqueza y disminuyeron en diversidad filogenética en niveles intermedios de perturbación forestal. Los endófitos de raíces no mostraron tal variación, pero ambas comunidades de endófitos difirieron en la composición de especies en función de la perturbación en el bosque. Las ubicaciones con menos perturbaciones generalmente tenían conjuntos distintivos de endófitos foliares, mientras que las ubicaciones más perturbadas tenían comunidades de endófitos regionalmente más homogéneas. Nuestros resultados muestran que los cambios en la riqueza de endófitos y la diversidad filogenética se pueden detectar en los tejidos de las hojas más de un año después de grandes tormentas. A su vez, los cambios generalizados en la composición de la comunidad de endófitos de hojas y raíces sugieren un efecto sutil y duradero de los huracanes que amerita un mayor estudio, lo que podría contribuir a la promoción de ensamblajes de comunidades de endófitos heterogéneas en diversos bosques. associated with plants may respond to forest disturbance. For example, soil microbes and fungi shift in community composition in response to disturbance by hurricanes, typically decreasing in diversity due to dominance by a smaller number of species (Cantrell et al., 2014;Eaton et al., 2020;Lodge, 1997). Soilborne microbial biomass also changes in response to forest damage, reflecting an increase in C and N concentrations in soil and altered soil processes (Lodge et al., 2016;Reed et al., 2020).
In comparison, relatively little is known about how forest disturbance due to hurricanes may impact endophyte communities.
Endophytes of forest plants typically are transmitted via air spora (leaves) and infection from soil (roots). Thus, the endophytes present in a tropical forest tree represent a subset of regionally available species, with strong impacts from ultraviolet radiation and other environmental factors that often are changed in disturbed forests (Oita et al., 2021). Because such endophytes typically have life stages outside of living tissues (see Herre et al., 2007;Nelson et al., 2020;Thomas et al., 2016), they may be affected by environmental shifts much like soil-inhabiting fungi. However, indirect effects of hurricanes on endophytes also may occur: hurricanes typically open forest canopies and increase litter accumulation, resulting in higher temperatures, changes in the light regime, and shifts in nutrient availability in the understory ) that impact plant physiology (Wen et al., 2008). Thus, forest disturbance due to hurricanes may select for: (1) endophyte species that are well adapted to withstand abiotic stress and/or a dynamic environment; (2) species that are more likely to have pathogenic or saprotrophic life stages, capitalizing on plant damage; and/or (3) specialized functional groups that facilitate the growth of host plants, ultimately benefiting the endophytes themselves (see Kandalepas, 2012;Maillard et al., 2020; see also Busby et al., 2022). If the patchy disturbance regimes imposed by hurricanes lead to different environmental pressures in different areas (see Kandalepas, 2012), such storms may lead to small-scale variation in communities at landscape scales (Gannon & Martin, 2014), ultimately shaping the distribution and subsequent impact of endophytes forest-wide.
The goal of this study was to understand how disturbance due to hurricanes may impact endophyte communities in forest plants.
Because impacts of hurricanes vary in intensity even within the same region (e.g., Lugo, 2000), it is possible to quantify effects of hurricanes on endophytes without confounding effects of large geographic distances or major shifts in edaphic or climate conditions.
Here, we measured disturbance due to Hurricane Maria in Carite State Forest of Puerto Rico ca. 16 months after the storm's passage, when permits and safe access were allowed by forest managers. In three sites, each comprised of three plots that represented a local gradient of disturbance due to the hurricane, we evaluated soil chemistry and collected leaves and roots, which we used in culture-free analyses to measure the richness, phylogenetic diversity, and composition of endophyte assemblages. We considered it plausible that these aspects of endophyte communities would vary linearly with the severity of forest damage, but also realized that relationships between endophyte richness or phylogenetic diversity and forest damage might not be linear, instead peaking or decreasing markedly only at intermediate levels of disturbance. Although often debated for its theoretical basis and empirical insights (e.g., see Fox, 2013asee Fox, , 2013bSheil & Burslem, 2013; see also Wang & Lin, 2019), the typical framework for considering such patterns is the intermediate disturbance hypothesis (IDH; Connell, 1978), which provides an expectation that richness or diversity is maximized when ecological disturbances are moderate.
Under this framework, we predicted that species richness of endophytes would be greatest where hurricane disturbances were intermediate (prediction 1). We also expected that endophytes in low-and high-disturbance sites might represent a broad phylogenetic diversity, with the former representing diverse, mature communities and the latter representing early successional communities.
Therefore, we predicted that phylogenetic diversity would decrease in intermediately disturbed areas relative to less-or more severely disturbed areas (prediction 2). Based on previous research showing that endophyte communities can vary spatially over small distances in neotropical forests (e.g., Higgins et al., 2014), we anticipated that endophyte communities would differ spatially but also as a function of forest disturbance (prediction 3). Finally, we predicted that areas with more extensive disturbance would have relatively similar endophyte communities to one another, whereas areas with less damage would have locally distinctive endophyte communities (prediction 4).
We selected three sites in humid and very humid subtropical forest ranging from 560 to 655 m a.s.l.: two areas at Charco Azul Recreational Area, and one at Guavate Recreational Area (Figure 1a).
The sites were similar in terms of forest structure and land use history prior to Hurricane Maria (H. Serrano, pers. comm.).
In each site, we characterized disturbance caused by Hurricane Maria following Basnet et al. (1992) with minor modifications. obtained for each category to obtain a forest damage score for each plot, which we log-transformed prior to analysis and considered representative of the relative disturbance in each site that could be ascribed to Hurricane Maria. Damage scores at the nine study plots ranged from 1 to 30 (Table 1).
In each plot we selected five representative plants. From each plant we collected 10 leaves and three root cores (6.3 cm diameter and ~10 cm depth, collected beneath the edge of the canopy of each individual). Root identity was confirmed by tracing roots where possible, and/or matching DNA sequences from leaves and roots as described below. We measured stem richness and density, canopy height, diameter at breast height (DBH; ≥1 cm) of the focal tree, and canopy cover in each plot as in Oita et al. (2021) (Table S1). For soil chemistry analyses, we collected three soil cores (6.3 cm diameter) from three randomly chosen corners of each plot after leaf litter was cleared. They were pooled by plot and sent to the Puerto Rico Tropical Agriculture Research Station (TARS) for quantification of nitrate, available phosphorus and potassium.

| Sample preparation
Soil samples were sieved to extract roots following Mukerji et al. (2002). Healthy leaves and roots collected from each individual plant were processed separately. In each case, tissues were cut into 1 cm-segments, rinsed in running tap water for 60 s, and surfacesterilized by agitation in 95% ethanol (30 s), 0.5% NaOCl − (2 min), and 70% ethanol (2 min; Arnold & Lutzoni, 2007). We randomly selected three segments from each sample and stored them in an Eppendorf tube containing 1 ml of cetyl trimethylammonium bromide (CTAB) buffer. This was repeated three times per individual plant, such that our final sampling included 9 cm 2 of leaf tissue and approximately 9 cm of root tissue per individual. Samples in CTAB were stored at −20°C prior to DNA extraction.

| DNA extraction, amplification, and sequencing for endophytes
We extracted total genomic DNA from leaves and roots with the Qiagen DNeasy Plant Mini Kit (Qiagen) following U'Ren and Arnold (2017a). Each CTAB tube represented a single DNA extraction, such that each individual was represented by three DNA extractions from leaf tissue, and three DNA extractions from root tissue. PCR products were stained with SYBR green and verified on a 2% agarose gel under ultraviolet (UV) light. The three PCR products per F I G U R E 1 (a) Carite State Forest, Puerto Rico. Each site included three plots representing a local gradient of forest damage scores, which were used to estimate disturbance due to Hurricane Maria. (b) Representative damage caused by Hurricane Maria at Carite State Forest. Left: plot with a low damage score, viewed from beyond the forest edge (site 1, plot 1). Right: plot with a high damage score (site 3, plot 3). Note: Plant taxa collected in each plot are shown in Table S4.

TA B L E 1
Abbreviations: m a.s.l., meters above sea level; ND, not determined.
sample were pooled and diluted with molecular grade water based on band intensity (U'Ren & Arnold, 2017a). Negative controls (water as template) and extraction blanks were treated as above for PCR1, and their products were processed in the second PCR (PCR2), in parallel with our samples.
In PCR2, we added adapters and unique barcodes to amplicons from PCR1 (see Oita et al., 2021;Sarmiento et al., 2017;U'Ren et al., 2019). Amplifications were carried out as above with a total volume of 20 μl (10 μl Table S2). The second mock community consisted of tiered DNA concentrations of the same fungal taxa (Daru et al., 2018; Table S2). Mock communities were treated as above and included with the samples for sequencing. We found that the even mock community produced unusual results consistent with contamination and degradation of DNA stocks, whereas the tiered mock community samples were effective positive controls. We therefore focused only on the tiered mock community as our positive control, as described below.

| Plant identification
Plants were identified based on Axelrod (2011), Little et al. (1974), the database of the Institute for Regional Conservation (https://www. to identify plants in cases in which morphological determination was not possible due to degradation of tissue following collection, when limited facilities were available for processing specimens. We amplified and sequenced the plant nuclear ribosomal ITS region because it was straightforward to amplify, but due to limitations in resolution we used ITS data to assign plants to the genus or family level only. PCR reaction volumes were 20 μl, including 10 μl of REDExtract-N-Amp PCR Reaction Mix (Sigma-Aldrich), 0.8 μl of each 10 μM primer (ITS5/ITS4; see White et al., 1990), 7.4 μl of molecular biology grade water (Fisher Scientific), and 1 μl of DNA template (Table S3). The cycling parameters were 3 min at 94°C, 36 cycles of 30 s denaturation at 94°C, 30 s annealing at 54°C, and 1 min extension at 72°C; and a final extension of 10 min at 72°C (Hoffman et al., 2008). PCR amplifications were verified by electrophoresis as described above.
Amplicons were treated with 1 μl ExoSAP-IT (Affymetrix) and incubated for 60 min at 37°C and 15 min at 80°C (Bell, 2008) prior to bidirectional Sanger sequencing at the University of Arizona Genetics Core. Sequencing reactions included 5 μM of the PCR primers and were processed on a 373xl DNA Sequencer (Applied Biosystems).
Contigs were assembled and edited with phred and phrap as implemented in Mesquite v.2.75 and verified by eye in Sequencher v.4.5 (Ewing & Green, 1998;Maddison & Maddison, 2017). Sequence data were compared via BLASTn to plant records in GenBank and crossreferenced with the herbarium and field guides mentioned above (Altschul et al., 1990; Table S4).
Samples that failed to be identified were amplified with alternate primer pairs that also target the plant nuclear ribosomal ITS region (Table S3)  Samples were amplified, sequenced, and identified as stated above.
DNA vouchers are retained at the University of Arizona.
Representative sequence data have been submitted to GenBank (accessions OP390837-OP390867).

| Illumina sequence editing and quality control
Illumina reads were demultiplexed at the IBEST Genomics Core.
Fast QC, Multi-QC, and fastq_eestats2 command in USEARCH were used to assess read quality and determine maxEE and cutoffs for forward and reverse reads Daru et al., 2018;Murray et al., 2015). Based on evaluations of the tiered mock communities, we selected forward reads for further evaluation, with a length cutoff of 240 bp and maxEE of 1.00.
Following these parameters, reads were trimmed and filtered with the -fastq_filter command in USEARCH (see . We used UNOISE3 and USEARCH to cluster sequences into Operational Taxonomic Units (OTUs) at 95% sequence similarity (UNOISE: -unoise3 command; USEARCH: -cluster_smallmem; see Daru et al., 2018). Reads from negative controls were treated as above and filtered from the remaining OTUs.
After filtering as above, the data set from the even mock community comprised 494,112 reads, and the data set from the tiered mock community comprised 390,333 reads. Following our initial quality control steps, 77 OTUs were detected in the mock community data sets when OTUs were based on 95% sequence similarity, yet only 31 were expected given the DNA inputs (Table S2). We found that most of the unexpected OTUs were in the even mock community samples, and that many species that should have been detected were not, leading us to conclude that the even mock community had been contaminated and was compromised. We therefore excluded the even mock community samples, instead examining the data for the tiered mock community. We found that by removing OTUs contributing <0.15% of the total read number (see Oita et al., 2021), we detected 25 OTUs among 355,185 reads for the tiered mock community.
With this basis, we implemented the same rules for the endophyte data sets (cutoff 240 bp; maxEE = 1.00, retaining OTUs based on 95% sequence similarity as long as they represented ≥0.15% of total read abundance). Ultimately this resulted in removal of OTUs with <13 occurrences (total reads). Leaf and root samples that yielded <100 reads also were removed prior to statistical analyses.
The final data sets consisted of 612,335 reads and 1094 OTUs representing endophytes from leaves, and 983,644 reads and 1128 OTUs representing endophytes from roots. Analyses of the tiered mock community revealed that read number and operon number were positively correlated ( Figure S1), such that we included abundance-based metrics in subsequent analyses.

| Statistical analyses
The purpose of our sampling design was to capture standing variation in leaf-and root endophyte communities as a function of forest disturbance due to Hurricane Maria. The heterogeneity of the forest required collection of different plant taxa in each plot. In turn, multiple plots per site, and multiple sites, were required to capture a range of damage levels across Carite State Forest.
Because plant species differed among study plots, our approach could confound damage score and plant identity as an explanatory factor in defining variation in endophyte communities. Therefore, we considered our sampling carefully prior to designing our statistical approaches. After the data processing and identification pro- We therefore concluded that there was not a systemic bias in our data due to the presence of different plant taxa on different plots.
This interpretation was upheld when we considered variation in endophyte richness among plant families: after variation in damage scores for plots was taken into account, there was no significant variation in endophyte richness among plant families for leaves or roots (p = .1382, p = .9215, respectively).  Figure S2).
With these provisions, we used multiple regression to consider how species richness of leaf and root endophytes differed as a func- Sordariomycetes, and Eurotiomycetes ( Figure S3), with indicator species referenced below.

| Endophyte richness in leaves, but not roots, peaked at intermediate disturbance
When analyzed independent of forest damage scores, the richness of endophytes in leaves and roots did not vary solely with plant family or PC soil ( Figures S4 and S5). We also found no meaningful variation in PC soil as a function of forest damage ( Figure S6).
However, when PC soil was taken into account and plant family was treated as a random factor, we found that endophyte richness in leaves varied as a function of forest damage scores (R 2 = .14, p = .0566), with a peak at intermediate disturbance levels (Table 2, Figure 2a). Endophyte richness in roots did not vary meaningfully as a function of forest damage scores (R 2 = .03, p = .6163; Table 2, Figure 2b). Thus data for foliar endophytes confirmed our first prediction, but endophytes belowground did not.  Figure 3a). However, a similar pattern was not detected for endophytes from roots (F 2,864 = 1.7120, p = .1811; Figure 3b).  Table S5. p = .2588) and roots (F 3,863 = 0.5438, p = .6524) did not vary among plots when variation due to forest damage was accounted for (data not shown).

| Endophyte communities varied spatially and as a function of forest disturbance
Sordariomycetes, Dothideomycetes, and Eurotiomycetes were the most common classes in all sites, and their prevalence was not sensitive to disturbance ( Figure S3). However, as anticipated (

| Endophyte communities in leaves, but not roots, became less distinctive in more-disturbed areas
As expected (prediction 4), plots with lower damage scores typically had leaf endophyte communities that were distinctive relative to those in other areas, whereas plots with more severe damage F I G U R E 5 Communities of root endophytes differed in species composition among sites and as a function of disturbance, as revealed by non-metric multidimensional scaling based on (a) Jaccard index and (b) Morisita index. Symbols show local disturbance: x, lowest damage score within a site; open circle, intermediate damage score within a site; closed circle, highest damage score within a site. Colors: black, site 1; aqua, site 2; blue, site 3. We excluded samples 135 and 132 (Jaccard) and 234 and 314 (Morisita; site, plot, and plant numbers as shown in Figure 4 and Table S4) because their endophytes were markedly different from all other samples. For statistical analyses see Table S5.
had leaf endophyte communities that were less distinctive relative to those in other sites (R 2 = .10; p < .0001; Table S7, Figure 6a).
The most common leaf-and root endophytes differed in their prevalence as a function of forest damage ( Figure S7). For leaf endophytes, 13 OTUs were identified as indicator species ( Figure S7). Of these, seven were identified in plots with intermediate disturbance

| DISCUSS ION
Hurricanes have profound effects on forest composition and ecosystem function in the Caribbean region (Tanner et al., 1991).
However, impacts of these disturbances on plant-microbe symbioses and related processes are not well understood. As hurricanes continue to increase in intensity and frequency in the Caribbean and western Atlantic (Kossin et al., 2017), cryptic effects on the symbionts with which forest trees interact may have long-term effects on forest regeneration and resilience (see Busby et al., 2022). Charting and their communities resembled the evolutionarily rich, mature endophyte assemblages of intact forest areas (see Arnold & Lutzoni, 2007).
In areas with high disturbance, leaf endophyte communities had high phylogenetic diversity but less were distinctive relative to other local collections, consistent with perturbation and/or early succession following the hurricane's direct effects. In sites with intermediate damage, foliar endophyte communities were species-rich but not especially phylogenetically diverse, suggesting that many closely related species may be assorting ecologically. These patterns were not observed in root endophytes, suggesting that their communities may react more slowly to disturbance, that all root endophyte communities regardless of the relative amount of disturbance had been perturbed in a relatively uniform fashion, or that spatial variation in root  Table S5). The quadratic fit based on disturbance alone is significant for leaf endophytes (R 2 = .58; p < .0001) but not root endophytes (R 2 = .10, p = .1152).
aspects of these variable associations is an important topic for future research, as it would link changes in composition of endophyte communities due to disturbance from hurricanes to the dynamics of plant communities in damaged and recovering forests.
More broadly, our study suggests that to understand impacts of hurricanes on forest processes mediated by endophytes, the appropriate sampling periods or scope may differ for tissues above and belowground. By the time of our study, leaves had generally flushed once or twice since the passage of Hurricane Maria, but the dynamics of root growth or death were not known. We anticipated that leaves would be sensitive to disturbance more immediately, and that roots could experience disturbance over longer time scales than that studied here (see Kandalepas, 2012). Even so, the differences we observed in community composition of both leaf and root endophytes are consistent with meaningful impacts by this hurricane on symbiont communities in both compartments. Our study thus provides a baseline for surveys that may detect additional changes in endophyte communities as this forest recovers from Hurricane Maria, and a basis for monitoring change following future disturbance events.
As our study was only retrospective, it is possible that variation in disturbance among our plots could reflect topography and related factors, and thus the differences we observed among endophytes could be signals of factors other than damage per se. To overcome this issue, we surveyed plots located in three sites, each encompassing a local damage gradient. Further, we accounted for differences in soil chemistry and addressed potential variation due to the presence of different plants among plots. Overall, we attribute the observed relationships of endophyte richness, phylogenetic diversity, and distinctiveness (leaves) and community composition (leaves and roots) to damage from the hurricane, but we encourage broader sampling to capture background variation in endophyte composition across the landscape. We note that the major groups of endophytes found in sites with low-and intermediate damage, including the most common genera, resemble those investigated by Lodge et al. (1996) in Manilkara bidentata of Puerto Rico, suggesting that our study captured representative endophyte communities for the region.
Increases in solar radiation and litter drying as a result of hurricane disturbance can result in the replacement of less tolerant fungal species by disturbance-adapted species and the occurrence of homogenous fungal communities, as observed following Hurricane Hugo in Puerto Rico (Lodge & Cantrell, 1995;Vandermeer et al., 2000). We suggest that these factors may partly explain alterations in the endophyte communities observed in our findings in response to hurricane disturbance. Nonetheless, Hubbell et al. (1999) demonstrated that disturbance events would not necessarily lead to the replacement of less-tolerant species by disturbance-adapted species, but instead by species that are simply abundant enough at the right place and time to establish opportunitistically. Hence, we cannot rule out that the spatial structure of fungal endophytes following a major disturbance event will depend on species abundance and recovery in the region, rather than the effect of stress alone.
Overall, our results show that endophytes are sensitive to damage due to major hurricanes, and that such sensitivity can be detected well after a year post-disturbance, particularly in foliar endophyte communities. We anticipate that widespread and repeated disturbance consistent with the increasing impact of hurricanes could lead to local extirpation of distinctive endophyte species that may have co-evolved with key components of island floras, and which may be important for the resilience of island ecosystems. Such effects may be relatively immediate in leaves and relatively slow in roots, a prediction that could be evaluated by resampling our plots in subsequent years. Overall, our work contributes to understanding how symbiont communities that can influence plant growth and resilience are sensitive to major abiotic disturbances, which is important for predicting their function in plants under an altered climate and developing strategies to restore severely damaged ecosystems in diverse regions of the world.

ACK N OWLED G M ENTS
The project was funded by the National Geographic Society

CO N FLI C T O F I NTE R E S T
No potential conflict of interest was identified by the authors.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available